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Abstract 

We present a new method devised to overcome the intrinsic difficulties associated 
to the numerical simulations of the Tsallis statistics. We use a standard Metropolis 
Monte Carlo algorithm at a fictitious temperature T', combined with a numerical 
integration method for the calculation of the entropy in order to evaluate the actual 
temperature T. We illustrate the method by applying it to the 2d-Ising model using 
a standard reweighting technique. 
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1 Introduction 



In 1988 C. Tsallis proposed a thermostatistics formalism based on a non- 
extensive entropy definition [1]. In the most recent formulation [2] of the Tsallis 
Statistics (TS) in the canonical ensemble, at fixed temperature T, observables 
Oq are obtained as averages of microscopic functions Of. 

w 

Oq = J20,P, (1) 
1=1 
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while the entropy is given by 



In these expressions, Pi are the escort probabihties[3] for configuration i — 
1, . . . ,W with energy sf. 

p. - ^ = [1 - (1 - <;)e./r']^ 
•"Er-^."Er.,[i-(i-<;te/r'l* 



and the additional rule that = whenever 1 — (1 — q)ei/T' < 0. Here, to 
simplify notation, one introduces the auxiliary parameter 

w w 
T'^{l-q)Y.^^P^ + nY.Pl'')-' (4) 

j=l i=l 



The probabilities Pj, the entropy 5"^ and the mean value defining the observable 
Oq depend, besides the temperature T, on a parameter g, which measures 
the degree of non-extensivity of the TS. It is not possible, in general, to 
solve the previous equations to give explicit expressions for the probabilities 
Pi as a function of q and T. An exception being the limit — > 1 in which 
one recovers the Boltzmann-Gibbs Statistics (BGS): Pi oc exp(— ej/T) and 
Si — — Z^i^i -Pi lii(-Pi). For systems with a small number W of configurations 
it is possible to solve equations (3,4) itcrativcly starting from an initial ansatz 
for Pj, e.g. the Boltzmann-Gibbs expression. However, this method is not 
useful for a system with a moderately large number of configurations. 

The fact that one can not give explicit expressions for the probabilities Pj 
has hampered the development of numerical methods to perform the usual 
Metropolis Monte Carlo or Molecular Dynamics simulations of the TS for 
interacting systems. Very recently, however, methods based on the numerical 
calculation of the number of configurations with a given energy have allowed 
the direct calculation of the necessary averages [4, 5]. In this paper, we introduce 
a new and more direct method, based on the standard Metropolis algorithm 
combined with a numerical integration, which can be used in many cases to 
perform the thermodynamic averages involved in the TS. In the next section 
we describe the method in some detail, and in section 3 we show the results 
of the application of the method to the Ising model as well as some of the 
difficulties encountered. 
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2 The Monte Carlo method 



As mentioned before, the main problem to perform a numerical simulation 
of a system described by the TS at fixed temperature T is that we do not 
have at hand the solution for the probabilities Pi, since the nonlinear equa- 
tions (3,4) have no explicit solution for q ^ 1. For q = 1 (BGS), it is 
Pi = -Z^^ exp(— £j/T), and one can use a variety of Monte Carlo techniques 
for the numerical calculation of the averages, Eq.(l). For example: in the 
Metropolis Monte Carlo algorithm[6], one generates a change in the config- 
uration i ^ j and the new configuration j is accepted with a probability 
mm{l, Pj/Pi) = min[l, exp( — (e^ — £j)/T)]. Notice that the partition function 
Z cancels out in the calculation of the acceptance probabilities. Unfortunately, 
since for g 7^ 1 the probabilities Pi are not known as a function of T, there is 
no trivial generalization of the Monte Carlo method to perform the averages 
(1) at fixed temperature T. The method we propose in this paper works in 
two steps: (i) we perform Monte Carlo simulations at a fixed value of the aux- 
iliary parameter T'; (ii) we then use Eq.(4) in order to determine the physical 
temperature T. We describe now in detail both steps. 

(i) To perform a Metropolis Monte Carlo simulation of the TS at a fixed 
"fictitious temperature" T', one proposes a change in the configuration i ^ j 
and accepts this change with probability min(l, Pj/Pj) [7]. Using Eq.(3), one 
notices that the normalizing factor cancels out: 



accep 



mm 



' T' - (1 - q)e, 
T - (1 - q)ei 



q - 

1-9 



(5) 



(it is also understood that the acceptance probability is zero if the configura- 
tion j is such that T' — (1 — q)ej < 0). By using this Monte Carlo algorithm, 
one generates a sequence of M representative configurations which arc dis- 
tributed according to the probability Eq.(3). The statistical averages, Eq.(l), 
are then approximated by sample averages Og — J2^=i Ok/M and the errors 
computed in the standard way [6]. 

(ii) To perform the T' ^ T transformation, we invert Eq.(4) using Eqs. (1,2): 

r - jl ^ q)U,{T') 

1 + (1 - q)Sg{T) ^ ' 



In this expression, the energy Uq{T') is obtained in the Monte Carlo simulation 
which is performed at fixed T'. In order to compute the entropy Sq{T') we make 
use of the thermodynamic relation 1/T = dSg/dUq[8] which, using Eq.(6), can 
be integrated between two equilibrium states characterized by values Tq and 
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T' of the parameter, to yield: 



1 ri + (l-g)5^(T') 



/ 



(7) 



l-q [l + {l-q)S,{U) 



T'-{l-q)U, 



The temperature T is finally given by: 



T'-{1- q)U,{T') 
l + il-q)S,iT^ 



T' 



exp {q 



T' - (1 - q)Ug 




In summary, one performs a Monte Carlo simulation using the Metropolis 
acceptance probability corresponding to states with a fixed value of T', starting 
from some initial value Tq to the desired valued of T'. In these simulations 
one computes, using the standard Monte Carlo procedure, the sample value 
of the energy Uq{T'). Finally, in order to obtain the physical temperature 
T one uses the previous expression Eq.(8) where the integral is computed 
numerically. The initial value Tq must be some limiting value in which the 
entropy Sq{TQ) is known. This depends on the problem, although obvious 
candidates are the very high or very low temperature configurations. One 
could think that a disadvantage of the method is that the integration in Eq. (8) 
requires a large number of simulations at different fictitious temperatures to be 
able to perform the T' — > T transformation accurately. However, we will show 
that the use of the reweighting techniques [9] reduces drastically the number 
of simulations. 



3 The 2-dimensional Ising Model 

To illustrate the method described in the previous section we present results 
for the short range Ising Model, defined by the following Hamiltonian: 



where each of the N — L x L spin variables Si can take the values ±1, and 
the sum X^<jj> runs over all nearest neighbor sites on a 2-dimensional lattice 
with periodic boundary conditions. 

The Metropolis Monte Carlo simulation is performed, as usual, by randomly 
choosing one of the spins. Si, and proposing a change in configuration in which 



) 



(9) 



<i,j> 
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Fig. 1. The points show the raw data obtained from a Metropohs Monte Carlo 
simulation using Eq.(5). The lines connecting the points have been obtained by 
using a reweighting extrapolation method, see Eq.(lO). We use here q = 0.8. 

the spin flips its value, —S^. This change is accepted with a probability 
given by Eq.(5) whereas, if rejected, the spin keeps its old value. 

In the Fig.(l) we show the raw Monte Carlo data for L — 4,10,20,30, in 
the case q = 0.8. In order to produce this plot, we have used the Metropolis 
algorithm at the points marked by symbols in the plot. The lines joining the 
points in the Fig. (1) have been obtained by a reweighting technique [9], which 
allows the calculation of mean values of observables for a T' different from the 
actual where the simulation was performed. In the case of the TS and the 
probabilities given by Eq.(2), the reweighting is based upon the exact relation: 



E.O(e)i/(£;T^)[i 



-il-q)e/r 
-{l-q)e/Ti 



l-(l-9)£/T' i 



(10) 



where H{e] T^) is the histogram of all the energy values generated in the Monte 
Carlo run at T^. Here 0{s) is the micro canonical mean value at energy e of 
the observable O. Notice that in Eq.(lO) the sums run over the energy levels 
e at variance with Eq.(l) where the sums run over all configurations i. 

Once we have obtained the averages as a function of T' as shown in the figure 
(1) for the energy, the value of the temperature T is obtained from Eq.(8) 
by performing the indicated integration from the chosen temperature Tq. In 
the Ising model, we have used the limit Tq = for which the representative 
configurations are the two ground states, Ug — 0, and the entropy is: Sq{TQ ~ 
0) = ^. 

/ l-q 
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Fig. 2. The shadowed area is the result of the integration needed in Eq.(8) to 
perform the T' — > T transformation. We use here AT = 10 x 10, and q = 0.8. 



In Fig. (2) we plot the function that has to be integrated in order to perform 
the T' — > T transformation. The main goal of this figure is to show that the 
function to integrate is, at least for these values of the parameters and system 
sizes, a smooth function. For the actual integration, we have used Simpson's 
3/8 rule. Finally, in Fig. (3) we plot the resulting internal energy Uq as a 
function of the actual temperature T, for several values of the parameter q 
and different system sizes. Wc note that for g < 1 one obtains a hysteresis- 
like loop that induces an ambiguity in the actual value of the energy. This is a 
generic feature of the TS, which can be resolved by applying a minimum free 
energy criterion[10] to choose the most stable solution. We have observed that 
for q > 1 and large system sizes (this is not the case in Fig. 3) the Monte Carlo 
simulations in some temperature ranges near the ferromagnetic transition take 
a very long time to equilibrate, thus preventing us from performing a very 
accurate measurement. We believe this is a generic feature of any dynamical 
updating scheme one could use to simulate the TS. 

Finally, we want to remark that Eq.(8) can be combined with any other sim- 
ulation technique which performs a sampling of the configuration space ac- 
cording to the probabihty Eq.(3). For instance, one could use the Molecular 
Dynamic methods using thermostats at T' [11] to study the dynamic behavior 
of systems with a large number of degrees of freedom. 
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Fig. 3. The internal energy Uq plotted as a function of the actual Temperature T for 
three values of the parameter q. The sizes used here are AT = 20 x 20 for 5 = 0.8, 1.0 
and AT = 4 X 4 for g = 1.2. 
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